﻿using System;
using System.Collections.Generic;
using System.Linq;
using BSU.GraduateWork.Data.Model;
using BSU.Math.Algebra;
using BSU.Math.Geometry;

namespace BSU.GraduateWork.BL
{
    public class GridFreeMethodFor2D
    {
        #region Public properties

        public Problem2D Problem { get; set; }

        public bool AutoGeneratePoints { get; set; }

        #endregion

        #region Private properties

        private List<Point3D> PrevU { get; set; }

        #endregion

        #region Constructors and Initializers

        public GridFreeMethodFor2D(Problem2D problem): this(problem, true)
        {
        }

        public GridFreeMethodFor2D(Problem2D problem, bool autoGeneratePoints)
        {
            AutoGeneratePoints = autoGeneratePoints;

            Problem = problem;
            Problem.Solution = new Solution2D();
            Problem.ExactSolution = new Solution2D();

            if (AutoGeneratePoints)
            {
                Problem.Solution.U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2);
            }
        }

        public GridFreeMethodFor2D(int n,
                              double x1, double x2,
                              double y1, double y2,
                              double g, double epsilon, double h, double alpha)
        {
            Problem = new Problem2D
            {
                InputParameters = new InputParameters2D
                {
                    N = n,

                    X1 = x1,
                    X2 = x2,
                    Y1 = y1,
                    Y2 = y2,

                    G = g,
                    Epsilon = epsilon,
                    H = h,
                    Alpha = alpha
                },
                Solution = new Solution2D
                {
                    IterationsCount = 0
                }
            };

            Problem.Solution.U = GenerateRegularPoints(Problem.InputParameters.N, Problem.InputParameters.X1, Problem.InputParameters.X2, Problem.InputParameters.Y1, Problem.InputParameters.Y2);
        }

        public static double F(double x, double y)
        {
            return 1;
        }

        #endregion

        #region Private methods

        private Matrix GetM(Point3D p, IList<Point3D> neighboringPoints)
        {
            int n = neighboringPoints.Count;

            var m = new Matrix(n + 1, 6);
            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];

                m[i, 0] = 1;
                m[i, 1] = point.X - p.X;
                m[i, 2] = point.Y - p.Y;
                m[i, 3] = 0.5 * (point.X - p.X) * (point.X - p.X);
                m[i, 4] = (point.X - p.X) * (point.Y - p.Y);
                m[i, 5] = 0.5 * (point.Y - p.Y) * (point.Y - p.Y);
            }
            m[n, 0] = Problem.Type == ProblemType.TypeThree ? -1 : 0;
            m[n, 1] = 0; m[n, 2] = 0; 
            m[n, 3] = 1; m[n, 4] = 0; m[n, 5] = 1; 

            return m;
        }

        private Matrix GetW(Point3D p, IList<Point3D> neighboringPoints)
        {
            int n = neighboringPoints.Count;
            var w = new Matrix(n + 1, n + 1);
            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];
                double norm = (new Vector(point) - new Vector(p)).Norm();

                if (norm / Problem.InputParameters.H <= 1)
                {
                    w[i, i] = System.Math.Exp(-Problem.InputParameters.Alpha * norm * norm / (Problem.InputParameters.H * Problem.InputParameters.H));
                }
                else
                {
                    w[i, i] = 0;
                }
            }
            w[n, n] = 1;
            return w;
        }

        private Vector GetB(Point3D p, IList<Point3D> neighboringPoints)
        {
            int n = neighboringPoints.Count;
            var b = new Vector(n + 1);

            for (int i = 0; i < n; i++)
            {
                var point = neighboringPoints[i];
                b[i] = GetPoint(point, neighboringPoints).Z;
            }
            b[n] = F(p.X, p.Y);
            return b;
        }

        public void Solve()
        {
            var startTime = DateTime.Now;

            Problem.Solution.IterationsCount = 0;
            do
            {
                PrevU = new List<Point3D>(Problem.InputParameters.N);
                foreach (Point3D p in Problem.Solution.U)
                {
                    PrevU.Add(new Point3D(p));
                }

                for (int i = 0; i < Problem.Solution.U.Count; i++)
                {
                    Point3D point = Problem.Solution.U[i];
                    if (!point.IsBoundary)
                    {
                        var neighboringPoints = GetNeighboringPoints(point);

                        var m = GetM(point, neighboringPoints);
                        var w = GetW(point, neighboringPoints);
                        var b = GetB(point, neighboringPoints);

                        var mt = m.Transpose();
                        var tmp = mt * w;

                        var a = (tmp * m).Inverse() * tmp * b;
                        point.Z = a[0];
                    }
                    else
                    {
                        point.Z = Problem.InputParameters.G;
                    }
                }
                Problem.Solution.IterationsCount++;
            } while (!IsIterationsStopped());

            var endTime = DateTime.Now;

            Problem.Solution.ExecutionTime = endTime - startTime;
        }

        private bool IsIterationsStopped()
        {
            double summ1 = 0;
            double summ2 = 0;

            for (int i = 0; i < Problem.Solution.U.Count; i++)
            {
                summ1 += System.Math.Abs(Problem.Solution.U[i].Z - PrevU[i].Z);
                summ2 += System.Math.Abs(Problem.Solution.U[i].Z);
            }
            return summ1 / summ2 < Problem.InputParameters.Epsilon;
        }

        private Point3D GetPoint(Point3D point, IEnumerable<Point3D> points)
        {
            var list = points
                .Where(p => p.X == point.X && p.Y == point.Y)
                .ToList();

            if (list.Count > 0)
            {
                return list[0];
            }
            return null;
        }

        private List<Point3D> GenerateRegularPoints(int count, double x1, double x2, double y1, double y2)
        {
            var points = new List<Point3D>(count);

            int boundaryPointsCount = (int)System.Math.Pow(count, 1.0 / 2);

            double hx = (x2 - x1) / boundaryPointsCount;
            double hy = (y2 - y1) / boundaryPointsCount;

            for (int i = 0; i <= boundaryPointsCount; i++)
            {
                for (int j = 0; j <= boundaryPointsCount; j++)
                {
                        points.Add(new Point3D(x1 + i * hx, y1 + j * hy, 0,
                                               i == 0 || j == 0 ||
                                               i == boundaryPointsCount || j == boundaryPointsCount));
                }
            }

            return points;
        }

        private IList<Point3D> GetNeighboringPoints(Point3D point)
        {
            var list = Problem.Solution.U
                .Where(p => GetDistance(point, p) <= Problem.InputParameters.H &&
                            (p.X != point.X && p.Y != point.Y))
                .ToList();

            return list;
        }

        private static double GetDistance(Point3D a, Point3D b)
        {
            return System.Math.Sqrt((a.X - b.X) * (a.X - b.X) +
                                    (a.Y - b.Y) * (a.Y - b.Y));
        }

        #endregion
    }
}
